home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / pc / FAQSYS18.ZIP / FAQS.DAT / EBERT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-08-02  |  39.3 KB  |  1,241 lines
  1. /*
  2.  *************************************************************************
  3.  *                               Ebert.c                                 *
  4.  *************************************************************************
  5.  * This file contains C code for procedures in my chapters of "Texturing *
  6.  * and Modeling: A Procedural Approach". This software may be used for   *
  7.  * non-profit use as long as the original author is credited.  If you    *
  8.  * find any bugs in this software, send email.                           *
  9.  *************************************************************************
  10.  * Author: David S. Ebert (ebert@cs.umbc.edu)                            *
  11.  * Copyright 1994 David S. Ebert                                         *
  12.  *************************************************************************
  13. */
  14.  
  15. /*
  16.  *************************************************************************
  17.  *                           WRITE_NOISE.C                               *
  18.  *  This procedure generates a noise function file for solid texturing   *
  19.  *                      by David S. Ebert                                *
  20.  *************************************************************************
  21.  * This noise file needs to be read in to your program for the rest of   *
  22.  * the procedures in this file to work.                                  *
  23.  *************************************************************************
  24.  */
  25.  
  26. #include <math.h>
  27. #include <stdio.h>
  28. #include "ebert.h"
  29. #include "macro.h"
  30.  
  31. #define SIZE           64  
  32. #define SIZE_1         65  
  33. #define POW_TABLE_SIZE 10000
  34. #define RAMP_SIZE      200
  35. #define OFFSET_SIZE    1024
  36. #define SWIRL_FRAMES   126 
  37. #define SWIRL          0.04986655               /* 2*PI/126  */
  38. #define SWIRL_AMOUNT   0.04986655               /* 2*PI/126  */
  39. #define DOWN           0.011904762
  40. #define SPEED          .012
  41. #define DOWN_AMOUNT    .0095
  42. #define DOWN_SMOKE     .015
  43. #define RAD1           .10
  44. #define RAD2           .14
  45. #define SWIRL_SMOKE    .062831853              /* 2*PI/100  */
  46. #define SWIRL_FRAMES_SMOKE 100                 
  47.  
  48. double drand48();
  49.  
  50. float  calc_noise();
  51. float  turbulence();
  52. rgb_td marble_color();
  53. double ease();
  54. void   transform_XYZ();
  55.  
  56. main()
  57. {
  58.   float   tmp,u,v;
  59.   long    i,j, k, ii,jj,kk;
  60.   float   noise[SIZE+1][SIZE+1][SIZE+1];
  61.   FILE    *noise_file;
  62.  
  63.  noise_file = fopen("noise.data","w");
  64.  
  65.  for (i=0; i<SIZE; i++)
  66.    for (j=0; j<SIZE; j++)
  67.      for (k=0; k<SIZE; k++)
  68.        {
  69.          noise[i][j][k] = (float)drand48();
  70.        }
  71.  
  72. /* This is a hack, but it works. Remember this is only done once.
  73.  */
  74.  
  75.  for (i=0; i<SIZE+1; i++)
  76.    for (j=0; j<SIZE+1; j++)
  77.      for (k=0; k<SIZE+1; k++)
  78.        {
  79.          ii = (i == SIZE)? 0:  i;
  80.          jj = (j == SIZE)? 0:  j;
  81.          kk = (k == SIZE)? 0:  k;
  82.          noise[i][j][k] = noise[ii][jj][kk];
  83.        }
  84.  
  85.  fwrite(noise,sizeof(float),(SIZE+1)*(SIZE+1)*(SIZE+1), noise_file);
  86.  fclose(noise_file);
  87. }
  88.  
  89. /*
  90.  ***********************************************************************
  91.  * The procedures below use the noise file written by the above        *
  92.  * write_noise program. Read the noise file into a variable called     *
  93.  * noise. These procedures also assume that the renderer has a global  *
  94.  * int variable frame_num  that contains the current frame number.     *
  95.  ***********************************************************************
  96.  */
  97.  
  98.  
  99. float   noise[SIZE+1][SIZE+1][SIZE+1];
  100. int     frame_num;
  101. /*
  102.  ***********************************************************************
  103.  *                            Calc_noise                               *
  104.  ***********************************************************************
  105.  * This is basically how the trilinear interpolation works. I lerp     *
  106.  * down the left front edge of the cube, then the right front edge of  *
  107.  * the cube(p_l, p_r). Then I lerp down the left back and right back   *
  108.  * edges of the  cube (p_l2, p_r2). Then I lerp across the front face  *
  109.  * between p_l  and p_r (p_face1). The I lerp across the back face     *
  110.  * between p_l2 and p_r2 (p_face2). Now I lerp along the line between  *
  111.  * p_face1 and  p_face2.                                               *
  112.  ***********************************************************************
  113.  */
  114.  
  115. float
  116. calc_noise(pnt)
  117.      xyz_td  pnt;
  118. {
  119.   float  t1;
  120.   float  p_l,p_l2,    /* value lerped down left side of face1 & face 2 */
  121.          p_r,p_r2,    /* value lerped down left side of face1 & face 2 */
  122.          p_face1,     /* value lerped across face 1 (x-y plane ceil of z) */
  123.          p_face2,     /* value lerped across face 2 (x-y plane floor of z) */
  124.          p_final;     /* value lerped through cube (in z)                  */
  125.  
  126.   extern float noise[SIZE_1][SIZE_1][SIZE_1];
  127.   float   tnoise;
  128.   register int      x, y, z,px,py,pz;
  129.   
  130.   px = (int)pnt.x;
  131.   py = (int)pnt.y;
  132.   pz = (int)pnt.z;
  133.   x = px &(SIZE-1); /* make sure the values are in the table           */
  134.   y = py &(SIZE-1); /* Effectively, replicates the table thought space */
  135.   z = pz &(SIZE-1);
  136.   
  137.   t1 = pnt.y - py;
  138.   p_l  = noise[x][y][z+1]+t1*(noise[x][y+1][z+1]-noise[x][y][z+1]);  
  139.   p_r  =noise[x+1][y][z+1]+t1*(noise[x+1][y+1][z+1]-noise[x+1][y][z+1]);
  140.   p_l2 = noise[x][y][z]+ t1*( noise[x][y+1][z] - noise[x][y][z]);     
  141.   p_r2 = noise[x+1][y][z]+ t1*(noise[x+1][y+1][z] - noise[x+1][y][z]);
  142.   
  143.   t1 = pnt.x - px; 
  144.   p_face1 = p_l + t1 * (p_r - p_l);
  145.   p_face2 = p_l2 + t1 * (p_r2 -p_l2);
  146.   
  147.   t1 = pnt.z - pz;
  148.   p_final =  p_face2 + t1*(p_face1 -p_face2);
  149.   
  150.   return(p_final);
  151. }
  152.  
  153.  
  154. /*
  155.  **********************************************************************
  156.  *                   TURBULENCE                                       *
  157.  **********************************************************************
  158.  */
  159.  
  160. float turbulence(pnt, pixel_size)
  161.     xyz_td   pnt;
  162.     float    pixel_size;
  163. {
  164.   float t, scale;
  165.  
  166.   t=0;
  167.   for(scale=1.0; scale >pixel_size; scale/=2.0)
  168.     {
  169.       pnt.x = pnt.x/scale; pnt.y = pnt.y/scale; pnt.z = pnt.z/scale;
  170.       t+=calc_noise(pnt)* scale ; 
  171.     }
  172.   return(t);
  173. }
  174.  
  175.  
  176. basic_gas(pnt,density,parms)
  177.      xyz_td  pnt; 
  178.      float   *density,*parms;
  179. {
  180.   float turb;
  181.   int   i;
  182.   static float pow_table[POW_TABLE_SIZE];
  183.   static int calcd=1;
  184.   
  185.   if(calcd)
  186.     { calcd=0;
  187.       for(i=POW_TABLE_SIZE-1; i>=0; i--)
  188.     pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
  189.                   parms[1]*2.0,(double)parms[2]);
  190.     }
  191.   turb =turbulence(pnt);
  192.   *density = pow_table[(int)(turb*(.5*(POW_TABLE_SIZE-1)))];
  193. }
  194.  
  195. /* 
  196.  *
  197.  *change to make veins of gas 
  198.  *
  199.  */
  200. /*       turb =(1.0 +sin(turbulence(pnt)*M_PI*5))*.5;*/
  201.  
  202.  
  203.  
  204.  
  205. steam_slab1(pnt, pnt_world, density,parms, vol)
  206.          xyz_td  pnt, pnt_world;
  207.          float   *density,*parms;
  208.          vol_td  vol;
  209. {
  210.   float        turb, dist_sq,density_max;
  211.   int          i, indx;
  212.   xyz_td       diff;
  213.   static float pow_table[POW_TABLE_SIZE], ramp[RAMP_SIZE], 
  214.                offset[OFFSET_SIZE];
  215.   static int   calcd=1;
  216.  
  217.       if(calcd)
  218.         { calcd=0;
  219.           for(i=POW_TABLE_SIZE-1; i>=0; i--)
  220.              pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
  221.                              parms[1]*2.0,(double)parms[2]);
  222.           make_tables(ramp);
  223.         }
  224.       turb =turbulence(pnt,.1);
  225.       *density = pow_table[(int)(turb*0.5*(POW_TABLE_SIZE-1))];
  226.  
  227.      /* determine distance from center of the slab ^2. */
  228.       XYZ_SUB(diff,vol.shape.center, pnt_world);
  229.       dist_sq = DOT_XYZ(diff,diff);
  230.       density_max = dist_sq*vol.shape.inv_rad_sq.y; 
  231.       indx = (int)((pnt.x+pnt.y+pnt.z)*100) & (OFFSET_SIZE -1);
  232.       density_max += parms[3]*offset[indx];
  233.  
  234.       if(density_max >= .25) /* ramp off if > 25% from center */
  235.         { i = (density_max -.25)*4/3*RAMP_SIZE; /* get table index 0:RAMP_SIZE-1 */
  236.           i=MIN(i,RAMP_SIZE-1);
  237.           density_max = ramp[i];
  238.           *density *=density_max;
  239.         }
  240.      }
  241.  
  242.      make_tables(ramp, offset)
  243.           float *ramp, *offset;
  244.      {
  245.        int   i;
  246.        float dist, result;
  247.        srand48(42);
  248.        for(i=0; i < OFFSET_SIZE; i++)
  249.         {
  250.          offset[i]=(float)drand48();
  251.         }
  252.       for(i = 0; i < RAMP_SIZE; i++)
  253.         { dist =i/(RAMP_SIZE -1.0);
  254.           ramp[i]=(cos(dist*M_PI) +1.0)/2.0;
  255.         }
  256.      }
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264. /* Addition needed to sleam_slab1 to get final steam rising from a teacup
  265.  */
  266.  
  267. /* dist = pnt_world.y - vol.shape.center.y;
  268.  if(dist > 0.0)
  269.    { dist = (dist +offset[indx]*.1)*vol.shape.inv_rad.y;
  270.      if(dist > .05)
  271.        { offset2 = (dist -.05)*1.111111;
  272.          offset2 = 1 - (exp(offset2)-1.0)/1.718282;
  273.          offset2 *=parms[1];
  274.          *density *= offset;
  275.        }
  276.      }
  277.      */
  278.  
  279.  
  280. /*
  281.  *********************************************************************
  282.  *                           Smoke_stream                            *
  283.  *********************************************************************
  284.  * parms[1] = Maximum density value - density scaling factor         *
  285.  * parms[2] = height for 0 density (end of ramping it off)           *
  286.  * parms[3] = height to start adding turbulence                      *
  287.  * parms[4] = height(length) for maximum turbulence;                 *
  288.  * parms[5] = height to start ramping density off                    *
  289.  * parms[6] = center.y                                               *
  290.  * parms[7] = speed for rising                                       *
  291.  * parms[8] = radius                                                 *
  292.  * parms[9] = max radius of swirling                                 *
  293.  *********************************************************************
  294.  */
  295.  
  296. smoke_stream(pnt,density,parms, pnt_world, vol)
  297.      xyz_td  pnt, pnt_world;
  298.      float   *density,*parms;
  299.      vol_td  *vol;
  300. {
  301.  float           dist_sq;
  302.  extern float    offset[OFFSET_SIZE];
  303.  xyz_td          diff;
  304.  xyz_td          hel_path, new_path, direction2, center;
  305.  double          ease(), turb_amount, theta_swirl, cos_theta, sin_theta;
  306.  static int      calcd=1;
  307.  static float    cos_theta2, sin_theta2;
  308.  static xyz_td   bottom;
  309.  static double   rad_sq, max_turb_length, radius, big_radius,
  310.                  st_d_ramp, d_ramp_length, end_d_ramp,
  311.                  inv_max_turb_length;
  312.  double          height, fast_turb, t_ease, path_turb, rad_sq2;
  313.  
  314.  if(calcd)
  315.    { bottom.x=0; bottom.z = 0;
  316.      bottom.y = parms[6];
  317.      radius   = parms[8];
  318.      big_radius = parms[9];
  319.      rad_sq   = radius*radius;
  320.      max_turb_length = parms[4];
  321.      inv_max_turb_length = 1/max_turb_length;
  322.      st_d_ramp = parms[5];
  323.      end_d_ramp = parms[2];
  324.      d_ramp_length = end_d_ramp - st_d_ramp;
  325.      theta_swirl = 45.0*M_PI/180.0; /* swirling effect */
  326.      cos_theta  = cos(theta_swirl);
  327.      sin_theta  = sin(theta_swirl);
  328.      cos_theta2 = .01*cos_theta;
  329.      sin_theta2 = .0075*sin_theta;
  330.      calcd=0;
  331.    }
  332.  
  333.  
  334.  height = pnt_world.y - bottom.y + calc_noise(pnt)*radius;
  335.    /* We don't want smoke below the bottom of the column */
  336.  if(height < 0)
  337.    { *density =0;  return;}
  338.  height -= parms[3];
  339.  if (height < 0.0)
  340.    height =0.0;
  341.  
  342.  /* calculate the eased turbulence, taking into account the value may be
  343.   *  greater than 1, which ease won't handle.
  344.   */
  345.  t_ease = height* inv_max_turb_length;
  346.  if(t_ease > 1.0)
  347.    { t_ease = ((int) (t_ease)) + ease( (t_ease - ((int)t_ease)), .001, .999);
  348.      if( t_ease > 2.5)
  349.        t_ease = 2.5;
  350.    }
  351.  else
  352.    t_ease = ease(t_ease, .5, .999); 
  353.  
  354.  /* Calculate the amount of turbulence to add in */    
  355.  fast_turb= turbulence(pnt, .1);
  356.  turb_amount = (fast_turb -0.875)* (.2 + .8*t_ease);
  357.  path_turb = fast_turb*(.2 + .8*t_ease);
  358.  
  359.  /* add turbulence to the height and see if it is above the top */
  360.  height +=0.1*turb_amount;
  361.  if(height > end_d_ramp)
  362.    { *density=0; return; }
  363.  
  364.  /* increase the radius of the column as the smoke rises */
  365.  if(height <=0)
  366.    rad_sq2 = rad_sq*.25;
  367.  else if (height <=end_d_ramp)
  368.    { rad_sq2 = (.5 + .5*(ease( height/(1.75*end_d_ramp), .5, .5)))*radius;
  369.      rad_sq2 *=rad_sq2;
  370.    }
  371.  
  372. /*
  373.  **************************************************************************
  374.  * move along a helical path 
  375.  **************************************************************************
  376.  *
  377.  
  378.  /* 
  379.   * calculate the path based on the unperturbed flow: helical path 
  380.   */
  381.  hel_path.x = cos_theta2 *(1+ path_turb)*(1+cos(pnt_world.y*M_PI*2)*.11)
  382.    *(1+ t_ease*.1)  + big_radius*path_turb;
  383.  hel_path.z = sin_theta2 *(1+path_turb)*
  384.    (1+sin(pnt_world.y*M_PI*2)*.085)* (1+ t_ease*.1) + .03*path_turb;
  385.  hel_path.y =  - path_turb;
  386.  XYZ_ADD(direction2, pnt_world, hel_path);
  387.  
  388.  /* adjusting the center point for ramping off the density based on the 
  389.   * turbulence of the moved point 
  390.   */
  391.  turb_amount *= big_radius;
  392.  center.x = bottom.x - turb_amount;
  393.  center.z = bottom.z + .75*turb_amount; 
  394.  
  395.  /* calculate the radial distance from the center and ramp off the density
  396.   * based on this distance squared.
  397.   */
  398.  diff.x = center.x - direction2.x;
  399.  diff.z = center.z - direction2.z;
  400.  dist_sq = diff.x*diff.x + diff.z*diff.z;
  401.  if(dist_sq > rad_sq2)
  402.    {*density=0; return;}
  403.  *density = (1-dist_sq/rad_sq2 + fast_turb*.05) * parms[1];
  404. if(height > st_d_ramp)
  405.    *density *=  (1- ease( (height - st_d_ramp)/(d_ramp_length), .5 , .5));
  406. }
  407.  
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  
  418.  
  419.      
  420. rgb_td marble(pnt)
  421.      xyz_td  pnt;
  422. {
  423.   float y;
  424.   y = pnt.y + 3.0*turbulence(pnt, .0125);
  425.   y = sin(y*M_PI); 
  426.   return (marble_color(y));
  427. }
  428.  
  429.  
  430. rgb_td marble_color(x)
  431.      float x;
  432. {
  433.   rgb_td  clr;
  434.   x = sqrt(x+1.0)*.7071; 
  435.   clr.g = .30 + .8*x;
  436.   x=sqrt(x);
  437.   clr.r = .30 + .6*x;
  438.   clr.b = .60 + .4*x; 
  439.   return (clr);
  440. }
  441.  
  442.  
  443.  
  444.  
  445. rgb_td marble_forming(pnt, frame_num, start_frame, end_frame)
  446.      xyz_td  pnt;
  447.      int     frame_num, start_frame, end_frame;
  448. {
  449.   float x, turb_percent, displacement;
  450.   
  451.   if(frame_num < start_frame)
  452.     { turb_percent=0;
  453.       displacement=0;
  454.     }
  455.   else if (frame_num >= end_frame)
  456.     { turb_percent=1;
  457.       displacement= 3;
  458.     }
  459.   else
  460.     { turb_percent= ((float)(frame_num-start_frame))/ (end_frame-start_frame);
  461.       displacement = 3*turb_percent;
  462.     }
  463.   
  464.   x = pnt.x + turb_percent*3.0*turbulence(pnt, .0125) - displacement;
  465.   x = sin(x*M_PI); 
  466.   return (marble_color(x));
  467. }
  468.  
  469. rgb_td marble_forming2(pnt, frame_num, start_frame, end_frame, heat_length)
  470.      xyz_td  pnt;
  471.      int     frame_num, start_frame, end_frame, heat_length;
  472. {
  473.   float       x, turb_percent, displacement, glow_percent;
  474.   rgb_td      m_color;
  475.   if(frame_num < (start_frame-heat_length/2) || 
  476.      frame_num > end_frame+heat_length/2)
  477.     glow_percent=0;
  478.   else if (frame_num < start_frame + heat_length/2)
  479.     glow_percent= 1.0 - ease( ((start_frame+heat_length/2-frame_num)/ 
  480.                    heat_length),0.4,0.6); 
  481.   else if (frame_num > end_frame-heat_length/2)
  482.     glow_percent =  ease( ((frame_num-(end_frame-heat_length/2))/ 
  483.                heat_length),0.4,0.6);
  484.   else 
  485.     glow_percent=1.0;
  486.   
  487.   if(frame_num < start_frame)
  488.     { turb_percent=0;
  489.       displacement=0;
  490.     }
  491.   else if (frame_num >= end_frame)
  492.     { turb_percent=1;
  493.       displacement= 3;
  494.     }
  495.   else
  496.     { turb_percent= ((float)(frame_num-start_frame))/(end_frame-start_frame);
  497.       turb_percent=ease(turb_percent, 0.3, 0.7);
  498.       displacement = 3*turb_percent;
  499.     }
  500.   
  501.   x = pnt.y + turb_percent*3.0*turbulence(pnt, .0125) - displacement;
  502.   x = sin(x*M_PI); 
  503.   m_color=marble_color(x);
  504.   glow_percent= .5* glow_percent;
  505.   m_color.r= glow_percent*(1.0)+ (1-glow_percent)*m_color.r;    
  506.   m_color.g= glow_percent*(0.4)+ (1-glow_percent)*m_color.g;    
  507.   m_color.b= glow_percent*(0.8)+ (1-glow_percent)*m_color.b;    
  508.   return(m_color);
  509. }
  510.  
  511. rgb_td moving_marble(pnt, frame_num)
  512.      xyz_td  pnt;
  513.      int     frame_num;
  514. {
  515.   float        x, tmp, tmp2;
  516.   static float down, theta, sin_theta, cos_theta;
  517.   xyz_td       hel_path, direction;
  518.   static int   calcd=1;
  519.   
  520.   if(calcd)
  521.     { theta =(frame_num%SWIRL_FRAMES)*SWIRL_AMOUNT; /* swirling effect */
  522.       cos_theta = RAD1 * cos(theta) + 0.5;
  523.       sin_theta = RAD2 * sin(theta) - 2.0;
  524.       down = (float)frame_num*DOWN_AMOUNT+2.0;
  525.       calcd=0;
  526.     }
  527.   tmp = calc_noise(pnt); /* add some randomness */
  528.   tmp2 = tmp*1.75;
  529.   
  530.   /* calculate the helical path */
  531.   hel_path.y = cos_theta + tmp;
  532.   hel_path.x = (- down)  + tmp2;
  533.   hel_path.z = sin_theta - tmp2;
  534.   XYZ_ADD(direction, pnt, hel_path);
  535.   
  536.   x = pnt.y + 3.0*turbulence(direction, .0125);
  537.   x = sin(x*M_PI); 
  538.   return (marble_color(x));
  539. }
  540.  
  541.  
  542.  
  543.  
  544. void fog(pnt, transp, frame_num)
  545.      xyz_td  pnt;
  546.      float   *transp;
  547.      int      frame_num;
  548. {
  549.   float tmp;
  550.   xyz_td direction,cyl;
  551.   double theta;
  552.   
  553.   pnt.x += 2.0 +turbulence(pnt, .1);
  554.   tmp = calc_noise(pnt);
  555.   pnt.y +=  4+tmp;
  556.   pnt.z += -2 - tmp; 
  557.   
  558.   theta =(frame_num%SWIRL_FRAMES)*SWIRL_AMOUNT; 
  559.   cyl.x =RAD1 * cos(theta); 
  560.   cyl.z =RAD2 * sin(theta);
  561.   
  562.   direction.x = pnt.x + cyl.x;
  563.   direction.y = pnt.y - frame_num*DOWN_AMOUNT;
  564.   direction.z = pnt.z + cyl.z;
  565.   
  566.   *transp = turbulence(direction, .015);
  567.   *transp = (1.0 -(*transp)*(*transp)*.275);
  568.   *transp =(*transp)*(*transp)*(*transp);
  569. }
  570.  
  571.  
  572. /*
  573.  ****************************************
  574.  * ANIMATING GASES
  575.  ****************************************
  576.  */
  577.  
  578. steam_moving(pnt, pnt_world, density,parms, vol)
  579.      xyz_td  pnt, pnt_world;
  580.      float   *density,*parms;
  581.      vol_td  vol;
  582. {
  583.   float  noise_amt,turb, dist_sq, density_max, offset2, theta, dist;
  584.   static float pow_table[POW_TABLE_SIZE], ramp[RAMP_SIZE], offset[OFFSET_SIZE];
  585.   extern int frame_num; 
  586.   xyz_td direction, diff;
  587.   int i, indx;
  588.   static int calcd=1;
  589.   static float down, cos_theta, sin_theta;
  590.  
  591.   if(calcd)
  592.     { calcd=0;
  593.       /* determine how to move the point through the space (helical path) */
  594.       theta =(frame_num%SWIRL_FRAMES)*SWIRL;       
  595.       down = (float)frame_num*DOWN*3.0 +4.0;
  596.       cos_theta = RAD1*cos(theta) +2.0;
  597.       sin_theta = RAD2*sin(theta) -2.0;
  598.       
  599.       for(i=POW_TABLE_SIZE-1; i>=0; i--)
  600.     pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
  601.                   parms[1]*2.0,(double)parms[2]);
  602.           make_tables(ramp);
  603.     }
  604.   
  605.   /* move the point along the helical path */
  606.   noise_amt = calc_noise(pnt);
  607.   direction.x = pnt.x + cos_theta + noise_amt;
  608.   direction.y = pnt.y - down + noise_amt;
  609.   direction.z = pnt.z +sin_theta + noise_amt;
  610.  
  611.   turb =turbulence(direction, .1);
  612.   *density = pow_table[(int)(turb*0.5*(POW_TABLE_SIZE-1))];
  613.   
  614.   /* determine distance from center of the slab ^2. */
  615.   XYZ_SUB(diff,vol.shape.center, pnt_world);
  616.   dist_sq = DOT_XYZ(diff,diff);
  617.   density_max = dist_sq*vol.shape.inv_rad_sq.y; 
  618.   indx = (int)((pnt.x+pnt.y+pnt.z)*100) & (OFFSET_SIZE -1);
  619.   density_max += parms[3]*offset[indx];
  620.   
  621.   if(density_max >= .25) /* ramp off if > 25% from center */
  622.     { i = (density_max -.25)*4/3*RAMP_SIZE; /* get table index 0:RAMP_SIZE-1 */
  623.       i=MIN(i,RAMP_SIZE-1);
  624.       density_max = ramp[i];
  625.       *density *=density_max;
  626.     }
  627.   
  628.   /* ramp it off vertically */
  629.   dist = pnt_world.y - vol.shape.center.y;
  630.   if(dist > 0.0)
  631.     { dist = (dist +offset[indx]*.1)*vol.shape.inv_rad.y;
  632.       if(dist > .05)
  633.     { offset2 = (dist -.05)*1.111111;
  634.       offset2 = 1 - (exp(offset2)-1.0)/1.718282;
  635.       offset2*=parms[1];
  636.       *density *= offset2;
  637.     }
  638.     }
  639. }
  640.  
  641.  
  642. volume_fog_animation(pnt, pnt_world, density, parms, vol)
  643.      xyz_td  pnt, pnt_world;
  644.      float   *density,*parms;
  645.      vol_td  *vol;
  646. {
  647.   float noise_amt, turb;
  648.    extern int frame_num; 
  649.    xyz_td direction;
  650.    int    indx;
  651.    static float pow_table[POW_TABLE_SIZE];
  652.    int    i;
  653.    static int calcd=1;
  654.    static float down, cos_theta, sin_theta, theta;
  655.  
  656.    if(calcd)
  657.      {
  658.        down = (float)frame_num*SPEED*1.5 +2.0;
  659.        theta =(frame_num%SWIRL_FRAMES)*SWIRL; /* get swirling effect */
  660.        cos_theta = cos(theta)*.1 + 0.5;       /* use a radius of .1  */
  661.        sin_theta = sin(theta)*.14 - 2.0;      /* use a radius of .14 */
  662.        calcd=0;
  663.        for(i=POW_TABLE_SIZE-1; i>=0; i--)
  664.          {
  665.           pow_table[i] = (float)pow(((double)(i))/(POW_TABLE_SIZE-1)*
  666.                          parms[1]*4.0,(double)parms[2]);
  667.          }
  668.      }
  669.      
  670.    /* make it move horizontally and add some noise to the movement  */
  671.    noise_amt = calc_noise(pnt); 
  672.    direction.x = pnt.x - down + noise_amt*1.5;
  673.    direction.y = pnt.y + cos_theta +noise_amt; 
  674.    direction.z = pnt.z + sin_theta -noise_amt*1.5;
  675.  
  676.    /* base the turbulence on the new point */
  677.    turb =turbulence(direction, .1);
  678.    *density = pow_table[(int)((turb*turb)*(.25*(POW_TABLE_SIZE-1)))];
  679.   
  680.    /* make sure density isn't greater than 1 */    
  681.    if(*density >1)
  682.      *density=1;
  683. }
  684.  
  685.  
  686.  
  687.  
  688. /*
  689.  *********************************************************************
  690.  *                           Rising_smoke_stream                     *
  691.  *********************************************************************
  692.  * parms[1] = Maximum density value - density scaling factor         *
  693.  * parms[2] = height for 0 density (end of ramping it off)           *
  694.  * parms[3] = height to start adding turbulence                      *
  695.  * parms[4] = height(length) for maximum turbulence;                 *
  696.  * parms[5] = height to start ramping density off                    *
  697.  * parms[6] = center.y                                               *
  698.  * parms[7] = speed for rising                                       *
  699.  * parms[8] = radius                                                 *
  700.  * parms[9] = max radius of swirling                                 *
  701.  *********************************************************************
  702.  */
  703.  
  704. rising_smoke_stream(pnt,density,parms, pnt_world, vol)
  705.      xyz_td  pnt, pnt_world;
  706.      float   *density,*parms;
  707.      vol_td  *vol;
  708. {
  709.   float  dist_sq;
  710.   extern float offset[OFFSET_SIZE];
  711.   extern int frame_num; 
  712.   static int calcd=1;
  713.   static float down, cos_theta2, sin_theta2;
  714.   xyz_td  hel_path, center, diff, direction2;
  715.   double ease(), turb_amount, theta_swirl, cos_theta, sin_theta;
  716.   static xyz_td   bottom;
  717.   static double   rad_sq, max_turb_length, radius, big_radius,
  718.                   st_d_ramp, d_ramp_length, end_d_ramp, inv_max_turb_length, 
  719.                   cos_theta3, down3, sin_theta3;
  720.   double          height, fast_turb, t_ease, path_turb, rad_sq2;
  721.  
  722.   if(calcd)
  723.     {
  724.       bottom.x = 0; bottom.z = 0;
  725.       bottom.y = parms[6];
  726.       radius   = parms[8];
  727.       big_radius = parms[9];
  728.       rad_sq   = radius*radius;
  729.       max_turb_length = parms[4];
  730.       inv_max_turb_length = 1/max_turb_length;
  731.       st_d_ramp = parms[5];
  732.       st_d_ramp =MIN(st_d_ramp, end_d_ramp);
  733.       end_d_ramp = parms[2];
  734.       d_ramp_length = end_d_ramp - st_d_ramp;
  735.       /* calculate the rotation about the helix axis based on the
  736.        * frame_number 
  737.        */
  738.     theta_swirl =(frame_num%SWIRL_FRAMES_SMOKE)*SWIRL_SMOKE;/*swirling effect*/
  739.     cos_theta  = cos(theta_swirl);
  740.     sin_theta  = sin(theta_swirl);
  741.     down = (float)(frame_num)*DOWN_SMOKE*.75 * parms[7];
  742.       /* Calculate sine and cosine of the different radii of the
  743.        * two helical helical paths
  744.        */
  745.     cos_theta2 = .01*cos_theta;
  746.     sin_theta2 = .0075*sin_theta;
  747.     cos_theta3= cos_theta2*2.25;
  748.     sin_theta3= sin_theta2*4.5;
  749.     down3= down*2.25;
  750.       calcd=0;
  751.     }
  752.  
  753.   height = pnt_world.y - bottom.y + calc_noise(pnt)*radius;
  754.     /* We don't want smoke below the bottom of the column */
  755.   if(height < 0)
  756.     { *density =0;  return;}
  757.   height -= parms[3];
  758.   if (height < 0.0)
  759.     height =0.0;
  760.  
  761.   /* calculate the eased turbulence, taking into account the value may be
  762.    *  greater than 1, which ease won't handle.
  763.    */
  764.   t_ease = height* inv_max_turb_length;
  765.   if(t_ease > 1.0)
  766.     { t_ease = ((int) (t_ease)) + ease( (t_ease - ((int)t_ease)), .001, .999);
  767.       if( t_ease > 2.5)
  768.         t_ease = 2.5;
  769.     }
  770.   else
  771.     t_ease = ease(t_ease, .5, .999); 
  772.  
  773.   /* move the point along the helical path before evaluating turbulence */
  774.  pnt.x += cos_theta3;
  775.  pnt.y -= down3;
  776.  pnt.z += sin_theta3;
  777.  
  778.   fast_turb= turbulence(pnt, .1);
  779.   turb_amount = (fast_turb -0.875)* (.2 + .8*t_ease);
  780.   path_turb = fast_turb*(.2 + .8*t_ease);
  781.  
  782.  /* add turbulence to the height and see if it is above the top */
  783.   height +=0.1*turb_amount;
  784.   if(height > end_d_ramp)
  785.    { *density=0; return; }
  786.  
  787.   /* increase the radius of the column as the smoke rises */
  788.   if(height <=0)
  789.     rad_sq2 = rad_sq*.25;
  790.   else if (height <=end_d_ramp)
  791.     {
  792.       rad_sq2 = (.5 + .5*(ease( height/(1.75*end_d_ramp), .5, .5)))*radius;
  793.       rad_sq2 *=rad_sq2;
  794.     }
  795.   else
  796.     rad_sq2 = rad_sq;
  797.  
  798.  /*
  799.   **************************************************************************
  800.   * move along a helical path plus add the ability to use tables
  801.   **************************************************************************
  802.   *
  803.  
  804.   /* 
  805.    * calculate the path base on the unperturbed flow: helical path 
  806.    */
  807.     hel_path.x = cos_theta2*(1+ path_turb)*(1+cos((pnt_world.y +down*.5)*
  808.         M_PI*2)    *.11)* (1+ t_ease*.1)  + big_radius*path_turb;
  809.   hel_path.z = sin_theta2 *(1+path_turb)*
  810.     (1+sin((pnt_world.y +down*.5)*M_PI*2)*.085)* (1+ t_ease*.1) + .03*path_turb;
  811.  
  812.  hel_path.y = (- down) - path_turb;
  813.   XYZ_ADD(direction2, pnt_world, hel_path);
  814.  
  815.   /* adjusting the center point for ramping off the density based on the 
  816.    *  turbulence of the moved point 
  817.    */
  818.   turb_amount *= big_radius;
  819.   center.x = bottom.x - turb_amount;
  820.   center.z = bottom.z + .75*turb_amount; 
  821.  
  822.  /* calculate the radial distance from the center and ramp off the density
  823.   * based on this distance squared.
  824.   */
  825.   diff.x = center.x - direction2.x;
  826.   diff.z = center.z - direction2.z;
  827.   dist_sq = diff.x*diff.x + diff.z*diff.z;
  828.  if(dist_sq > rad_sq2)
  829.    {*density=0; return;}
  830.  *density = (1-dist_sq/rad_sq2 + fast_turb*.05) * parms[1];
  831. if(height > st_d_ramp)
  832.    *density *=  (1- ease( (height - st_d_ramp)/(d_ramp_length), .5 , .5));
  833. }
  834.  
  835.  
  836.  
  837.  
  838.  
  839.  
  840.  
  841.  
  842.  
  843.  
  844.  
  845. spherical_attractor(point, ff, direction, density_scaling, velocity, 
  846.             percent_to_use)
  847.      xyz_td       point, *direction;
  848.      flow_func_td ff;
  849.      float        *density_scaling, *velocity, *percent_to_use;
  850. {
  851.   float        dist, d2;
  852.   
  853.   /*calculate distance and direction from center of attractor  */
  854.   XYZ_SUB(*direction, point, ff.center);
  855.   dist=sqrt(DOT_XYZ(*direction,*direction));
  856.   
  857.   /* set the density scaling and the velocity to 1 */
  858.   *density_scaling=1.0;
  859.   *velocity=1.0;
  860.   
  861.   /* calculate the falloff factor (cosine) */
  862.   if(dist > ff.distance)
  863.     *percent_to_use=0;
  864.   else if (dist < ff.falloff_start)
  865.     *percent_to_use=1.0;
  866.   else
  867.     { d2 = (dist - ff.falloff_start)/(ff.distance - ff.falloff_start);
  868.       *percent_to_use = (cos(d2*M_PI)+1.0)*.5;
  869.     }
  870. }
  871.  
  872. calc_vortex(pt, ff, direction,percent_to_use, frame_num, velocity)
  873.      xyz_td       *pt, *direction;
  874.      flow_func_td *ff;
  875.      float        *percent_to_use, *velocity;
  876.      int          frame_num;
  877. {
  878.   static tran_mat_td  mat={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
  879.        xyz_td              dir, pt2, diff;
  880.        float               theta, dist, d2, dist2;
  881.        float               cos_theta,sin_theta,compl_cos, ratio_mult;
  882.  
  883.         /*calculate distance from center of vortex */
  884.        XYZ_SUB(diff,(*pt), ff->center);
  885.        dist=sqrt(DOT_XYZ(diff,diff));
  886.        dist2 = dist/ff->distance;
  887.         /* calculate angle of rotation about the axis */
  888.        theta = (ff->parms[0]*(1+.001*(frame_num)))/
  889.                (pow((.1+dist2*.9), ff->parms[1])); 
  890.  
  891.         /* calculate the matrix for rotating about the cylinder's axis */
  892.        calc_rot_mat(theta, ff->axis, mat);
  893.        transform_XYZ((long)1,mat,pt,&pt2);
  894.        XYZ_SUB(dir,pt2,(*pt));
  895.        direction->x = dir.x;
  896.        direction->y = dir.y;
  897.        direction->z = dir.z;
  898.   
  899.         /* Have the maximum strength increase from frame parms[4] to 
  900.          * parms[5]to a maximum of parms[2]   */
  901.        if(frame_num < ff->parms[4])
  902.          ratio_mult=0;
  903.        else if (frame_num <= ff->parms[5])
  904.           ratio_mult = (frame_num - ff->parms[4])/
  905.                        (ff->parms[5] - ff->parms[4])* ff->parms[2];
  906.        else
  907.          ratio_mult = ff->parms[2];
  908.  
  909.         /* calculate the falloff factor */
  910.        if(dist > ff->distance)
  911.          { *percent_to_use=0;
  912.            *velocity=1;
  913.          }
  914.        else if (dist < ff->falloff_start)
  915.          { *percent_to_use=1.0 *ratio_mult;
  916.            /*calc velocity  */
  917.            *velocity= 1.0+(1.0 - (dist/ff->falloff_start));
  918.          }
  919.        else
  920.          { d2 = (dist - ff->falloff_start)/(ff->distance - ff->falloff_start);
  921.            *percent_to_use = (cos(d2*M_PI)+1.0)*.5*ratio_mult;
  922.            *velocity= 1.0+(1.0 - (dist/ff->falloff_start));
  923.          }
  924. }
  925.  
  926.  
  927.  
  928. /*
  929.  ***************************************************************************
  930.  * Move the Volume of the Steam 
  931.  ***************************************************************************
  932.  * the shifting is based on the height above the cup (parms[6]->parms[7])
  933.  * and the frame range for increasing the strength of the attractor.
  934.  *  This is gotten from ratio_mult that is calculated above in calc_vortex. 
  935.  ***************************************************************************
  936.  */
  937.  
  938. /* Have the maximum strength increase from frame parms[4] to 
  939.  * parms[5]to a maximum of parms[2]   
  940.  */
  941. /* if(frame_num < ff->parms[4])
  942.     ratio_mult=0;
  943.    else if (frame_num <= ff->parms[5])
  944.      ratio_mult = (frame_num - ff->parms[4])/(ff->parms[5] - ff->parms[4])
  945.      * ff->parms[2];
  946.      if(point.y < ff->parms[6])
  947.        x_disp=0;
  948.      else
  949.       { if(point.y <= ff->parms[7])
  950.           d2 =COS_ERP((point.y - ff->parms[6])/(ff->parms[7] -ff->parms[6]));
  951.         else 
  952.           d2=0;
  953.         x_disp = (1-d2)*ratio_mult*parms[8]+calc_noise(point)*ff->parms[9];
  954.       }
  955.     return(x_disp);
  956. */
  957.  
  958.  
  959.  
  960.  
  961.  
  962.  
  963.  
  964.  
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973.  
  974.  
  975. /*
  976.  *********************************************************************
  977.  * parms[1] = Maximum density value: density scaling factor          *
  978.  * parms[2] = exponent for  density scaling                          *
  979.  * parms[3] = x resolution for Perlin's trick (0-640)                *
  980.  * parms[8] = 1/radius of fuzzy area for perlin's trick (> 1.0)      *
  981.  *********************************************************************
  982.  */
  983.  
  984. molten_marble(pnt, density, parms,vol)
  985.      xyz_td  pnt;
  986.      float   *density,*parms;
  987.      vol_td  vol;
  988. {
  989.   float  parms_scalar, turb_amount;
  990.   
  991.   turb_amount = solid_txt(pnt,vol);
  992.   *density = (pow(turb_amount, parms[2]) )*0.35 +.65;
  993.   /* Introduce a harder surface quicker. parms[3] is multiplied by 1/640 */
  994.   *density *=parms[1]; 
  995.   parms_scalar = (parms[3]*.0015625)*parms[8];
  996.   *density= (*density-.5)*parms_scalar +.5;
  997.   *density = MAX(0.2, MIN(1.0,*density));
  998. }
  999.  
  1000. /* 
  1001.  *=====================================================================*
  1002.  * --------------------------------------------------------------------*
  1003.  * ease-in/ease-out                                                    *
  1004.  * --------------------------------------------------------------------*
  1005.  * By Dr. Richard E. Parent, The Ohio State University                 *
  1006.  * (parent@cis.ohio-state.edu)                                         *
  1007.  * --------------------------------------------------------------------*
  1008.  * using parabolic blending at the end points                          *
  1009.  * first leg has constant acceleration from 0 to v during time 0 to t1 *
  1010.  * second leg has constant velocity of v during time t1 to t2          *
  1011.  * third leg has constant deceleration from v to 0 during time t2 to 1 *
  1012.  * these are integrated to get the 'distance' traveled at any time     *
  1013.  * --------------------------------------------------------------------*
  1014.  */
  1015. double ease(t,t1,t2)
  1016. double  t,t1,t2;
  1017. {
  1018.   double  ans,s,a,b,c,nt,rt;
  1019.   double  v,a1,a2;
  1020.  
  1021.   v = 2/(1+t2-t1);  /* constant velocity attained */
  1022.   a1 = v/t1;        /* acceleration of first leg */
  1023.   a2 = -v/(1-t2);   /* deceleration of last leg */
  1024.  
  1025.   if (t<t1) {
  1026.     rt = 0.5*a1*t*t;       /* pos = 1/2 * acc * t*t */
  1027.   }
  1028.   else if (t<t2) {
  1029.     a = 0.5*a1*t1*t1;      /* distance from first leg */
  1030.     b = v*(t-t1);            /* distance = vel * time  of second leg */
  1031.     rt = a + b;
  1032.   }
  1033.   else {
  1034.     a = 0.5*a1*t1*t1;      /* distance from first leg */
  1035.     b = v*(t2-t1);           /* distance from second leg */
  1036.     c = ((v + v + (t-t2)*a2)/2) * (t-t2);  /* distance = ave vel. * time */
  1037.     rt = a + b + c;
  1038.   }
  1039.   return(rt);
  1040. }
  1041.  
  1042.  
  1043.  
  1044. /* transform_xyz: transform an array of xyz_td and output cooresponding
  1045.  * xyz_td pnts which are transformed by the matrix input. 
  1046.  * It should actually only be be used for scale, rotate, and translate.
  1047.  */
  1048.  
  1049. void
  1050. transform_XYZ(num_pts, mat, inpts, outpts)
  1051.  
  1052.     long        num_pts;
  1053.     tran_mat_td mat;
  1054.     xyz_td      *inpts;
  1055.     xyz_td     *outpts;
  1056.     
  1057. {
  1058.     long i;
  1059.     float w;
  1060.     xyz_td *tmp, tmp2;
  1061.  
  1062.     if(num_pts==1)
  1063.       {
  1064.     tmp2.x = mat[0][0]*inpts[0].x + mat[1][0]*inpts[0].y
  1065.              + mat[2][0]*inpts[0].z + mat[3][0];
  1066.     tmp2.y = mat[0][1]*inpts[0].x + mat[1][1]*inpts[0].y
  1067.                  + mat[2][1]*inpts[0].z + mat[3][1];
  1068.     tmp2.z = mat[0][2]*inpts[0].x + mat[1][2]*inpts[0].y
  1069.                  + mat[2][2]*inpts[0].z + mat[3][2];
  1070.     *outpts=tmp2;
  1071.     return;
  1072.       }
  1073.     /* else */
  1074.     tmp =(xyz_td *)malloc(sizeof(xyz_td)*num_pts);
  1075.  
  1076.     for (i = 0; i < num_pts; i++)
  1077.       {
  1078.     tmp[i].x = mat[0][0]*inpts[i].x + mat[1][0]*inpts[i].y
  1079.                 + mat[2][0]*inpts[i].z + mat[3][0];
  1080.     tmp[i].y = mat[0][1]*inpts[i].x + mat[1][1]*inpts[i].y
  1081.                 + mat[2][1]*inpts[i].z + mat[3][1];
  1082.     tmp[i].z = mat[0][2]*inpts[i].x + mat[1][2]*inpts[i].y
  1083.                     + mat[2][2]*inpts[i].z + mat[3][2];
  1084.     outpts[i].x = tmp[i].x;
  1085.     outpts[i].y = tmp[i].y;
  1086.     outpts[i].z = tmp[i].z;
  1087.       }
  1088.     
  1089. }
  1090.  
  1091. //************************ EBERT.H ****************************
  1092.  
  1093. /* EDGE standard header file for graphics programs
  1094.  * Stripped down version 3/11/94 
  1095.  */
  1096.  
  1097. #ifndef GRAPHICS_H
  1098. #define GRAPHICS_H
  1099. #endif
  1100.  
  1101. #define       ABS(a)          ((a) > 0 ? (a) : -(a))
  1102. double fabs(), sqrt(), sin(), cos(), tan();
  1103.  
  1104.   typedef struct xyz_td
  1105.           {
  1106.               float    x,
  1107.                        y,
  1108.                        z;
  1109.       } xyz_td;
  1110.  
  1111.  
  1112.  
  1113.   typedef struct rgb_td
  1114.           {
  1115.               float    r,
  1116.                        g,
  1117.                        b;
  1118.       } rgb_td;
  1119.  
  1120.  
  1121. typedef float tran_mat_td[4][4];
  1122.  
  1123.  
  1124. /*
  1125.  *****************************************************************************
  1126.  *                  FUNCTION FIELD DEFINITIONS
  1127.  *****************************************************************************
  1128.  */
  1129.  
  1130. typedef unsigned char flow_func_type;
  1131.  
  1132. typedef struct flow_func_td       
  1133.     {
  1134.       short         func_type;
  1135.       xyz_td        center, axis;
  1136.       float         distance;
  1137.       float         falloff_start;
  1138.       short         falloff_type;
  1139.       float         parms[20];
  1140.     } flow_func_td;
  1141.  
  1142. typedef struct vol_shape_td
  1143.     {
  1144.       xyz_td         center;
  1145.       xyz_td         inv_rad_sq;
  1146.       xyz_td         inv_rad;
  1147.     } vol_shape_td;
  1148.  
  1149. typedef struct vol_td
  1150.     {
  1151.       short          obj_type;
  1152.       float          y_obj_min,    /* min & max y per volume             */
  1153.                  y_obj_max;
  1154.       float          x_obj_min,    /* min & max x per volume             */
  1155.                  x_obj_max;
  1156.       int            shape_type;   /* 0=sphere, 1=box;                   */
  1157.       int            self_shadow;  /* 0=no, 1=yes                        */
  1158.       xyz_td         scale_obj;   /*how to scale obj to get into -1:1 space*/
  1159.       xyz_td         trans_obj;   /*how to trans obj to get into -1:1 space*/
  1160.       vol_shape_td   shape;
  1161.       int            funct_name;
  1162.       rgb_td         color;
  1163.       float          amb_coef;
  1164.       float          diff_coef;
  1165.       float          spec_coef;
  1166.       float          spec_power;
  1167.       int            illum_type;  /* illumination - phong, blinn, c-t, gas */
  1168.       int            color_type;  /* constant, solid*/
  1169.       float          indx_refrac;
  1170.     } vol_td;
  1171.  
  1172. //************************* MACRO.H **********************************
  1173.  
  1174. #define XYZ_ADD(c,a,b)   { (c).x = (a).x + (b).x; (c).y = (a).y + (b).y;\
  1175.                            (c).z = (a).z + (b).z;}
  1176. #define XYZ_SUB(c,a,b)   { (c).x = (a).x - (b).x; (c).y = (a).y - (b).y; \
  1177.                            (c).z = (a).z - (b).z;}
  1178. #define XYZ_MULT(c,a,b)  { (c).x = (a).x * (b).x; (c).y = (a).y * (b).y; \
  1179.                            (c).z = (a).z * (b).z;}
  1180. #define XYZ_DIV(c,a,b)   { (c).x = (a).x / (b).x; (c).y = (a).y / (b).y; \
  1181.                            (c).z = (a).z / (b).z;}
  1182.  
  1183. #define XYZ_INC(c,a)     { (c).x += (a).x; (c).y += (a).y;\
  1184.                            (c).z += (a).z;}
  1185. #define XYZ_DEC(c,a)     { (c).x -= (a).x; (c).y -= (a).y;\
  1186.                            (c).z -= (a).z;}
  1187. #define XYZ_ADDC(c,a,b)  { (c).x = (a).x + (b); (c).y = (a).y + (b);\
  1188.                            (c).z = (a).z + (b);}
  1189. #define XYZ_SUBC(c,a,b)  { (c).x = (a).x - (b); (c).y = (a).y - (b); \
  1190.                            (c).z = (a).z - (b);}
  1191. #define XYZ_MULTC(c,a,b) { (c).x = (a).x * (b); (c).y = (a).y * (b);\
  1192.                            (c).z = (a).z * (b);}
  1193. #define XYZ_DIVC(c,a,b)  { (c).x = (a).x / (b); (c).y = (a).y / (b);\
  1194.                            (c).z = (a).z / (b);}
  1195. #define XYZ_COPY(b,a)    { (b).x = (a).x; (b).y = (a).y; (b).z = (a).z; }
  1196. #define XYZ_COPYC(b,a)   { (b).x = (a); (b).y = (a); (b).z = (a); }
  1197. #define  DOT_XYZ(p1, p2)  ((p1).x * (p2).x + (p1).y * (p2).y + (p1).z * (p2).z)
  1198.  
  1199. #define CROSS_XYZ(c,a,b) { (c).x = (a).y * (b).z - (a).z * (b).y; \
  1200.                (c).y = (a).z * (b).x - (a).x * (b).z; \
  1201.                (c).z = (a).x * (b).y - (a).y * (b).x; }
  1202. #define CROSS_3(c,a,b)   { (c)[0] = (a)[1] * (b)[2] - (a)[2] * (b)[1]; \
  1203.                (c)[1] = (a)[2] * (b)[0] - (a)[0] * (b)[2]; \
  1204.                (c)[2] = (a)[0] * (b)[1] - (a)[1] * (b)[0]; }
  1205.  
  1206. #define NORMALIZE_XYZ(v)  { float __tmpnormval; \
  1207.                 __tmpnormval = (double)sqrt((v).x*(v).x +  \
  1208.                         (v).y*(v).y + (v).z*(v).z); \
  1209.                   (v).x /= __tmpnormval; \
  1210.                   (v).y /= __tmpnormval;  \
  1211.                   (v).z /= __tmpnormval;}
  1212. #define R_NORMALIZE_XYZ(v)  { float __tmpnormval; \
  1213.                 (((__tmpnormval = (double)sqrt((v).x*(v).x +  \
  1214.                         (v).y*(v).y + (v).z*(v).z)) \
  1215.                   == 0.0) ? FALSE : ((v).x /= __tmpnormval,  \
  1216.                          (v).y /= __tmpnormval,  \
  1217.                          (v).z /= __tmpnormval,TRUE));}
  1218.  
  1219. #define RGB_ADD(c,a,e)   { (c).r = (a).r + (e).r; (c).g = (a).g + (e).g; \
  1220.                            (c).b = (a).b + (e).b; }
  1221. #define RGB_SUB(c,a,e)   { (c).r = (a).r - (e).r; (c).g = (a).g - (e).g; \
  1222.                            (c).b = (a).b - (e).b; }
  1223. #define RGB_MULT(c,a,e)  { (c).r = (a).r * (e).r; (c).g = (a).g * (e).g; \
  1224.                            (c).b = (a).b * (e).b; }
  1225. #define RGB_ADDC(c,a,e)  { (c).r = (a).r + (e); (c).g = (a).g + (e); \
  1226.                            (c).b = (a).b + (e); }
  1227. #define RGB_SUBC(c,a,e)  { (c).r = (a).r - (e); (c).g = (a).g - (e); \
  1228.                            (c).b = (a).b - (e); }
  1229. #define RGB_MULTC(c,a,e) { (c).r = (a).r * (e); (c).g = (a).g * (e); \
  1230.                            (c).b = (a).b * (e); }
  1231. #define RGB_DIVC(c,a,e)  { (c).r = (a).r / (e); (c).g = (a).g / (e); \
  1232.                            (c).b = (a).b / (e); }
  1233. #define RGB_COPY(c,a)    { (c).r = (a).r; (c).g = (a).g; (c).b = (a).b; }
  1234. #define RGB_COPYC(c,a)   { (c).r = (a); (c).g = (a); (c).b = (a); }
  1235.  
  1236. #undef MIN
  1237. #undef MAX
  1238. #define    MIN(a, b)       ((a) < (b) ? (a) : (b))
  1239. #define    MAX(a, b)       ((a) > (b) ? (a) : (b))
  1240.  
  1241.